; docformat = 'rst'
;
;+
;
; :Purpose:
;   Calculate NIterations random realizations of vector x given covariance matrix P
;   
; :Inputs:
;   x: vector of mean parameters
;   P: covariance matrix of vector x
;   nIterations: number of random realisations of x
;
; :Keywords:
;
; :Outputs:
;   An n_elements(x) by nIterations array of random realisations of x
;   
; :Example::
;   xr=mr_cholesky_covariance(x, P, 100)
;   
; :History:
; 	Written by: Matt Rigby, University of Bristol, Jan 29, 2013
;
;-
function mr_cholesky_covariance, x, P, nIterations

  compile_opt idl2

  seed=systime(/seconds)
  nX=n_elements(x)
  xr=fltarr(nX, nIterations)
  L=double(P)
  choldc, l, ld, /double

  for n=0, n_elements(l[0, *])-1 do l[n:-1, n]=0.
  l+=diag_matrix(ld)
  for n=0, nIterations-1 do xr[*, n]=L##randomn(seed, nX) + x

  return, xr
  
end